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EXECUTIVE  SUMMARY 


A  new  formulation  for  the  calculation  of  radar  cross 
section  of  objects  of  arbitrary  shape  based  on  the  solution 
of  the  electromagnetic  wave  equation  in  a  natural  coordinate 
system  that  conforms  to  the  shape  of  the  boundaries  has  been 
developed.  The  natural  coordinate  system  is  determined  from 
the  solution  of  an  elliptical  partial  differential  equation 
subject  to  Dirichlet  boundary  conditions.  The  basic  approach 
is  to  solve  the  electromagnetic  wave  equation  expressed  in 
the  natural  coordinate  system  in  a  region  between  the 
scattering  object  and  an  fictitious  interface  that  fully 
envelopes  the  object  (a  circle  for  two  dimensional  objects 
and  a  sphere  for  three  dimensional  objects) .  Outside  the 
interface,  the  wave  field  can  be  expressed  as  infinite  series 
of  Bessel  types  of  functions.  A  unique  solution  is  obtained 
by  matching  the  boundary  conditions  at  the  fictitious 
interface.  -  ,  ^  ^ 

Under  the  six-months  DOD/SBIR  Phase  I  Prograun  funded  by 
DARPA  ,  the  development  of  an  algorithm  based  on  this 
formulation  to  calculate  the  radar  cross  section  of  two 
dimensional  objects  to  demonstrate  the  advantages  of  this 
approach  has  been  successfully  carried  out.  The  formulation 
is  completely  general  and  can  be  extended  to  analyze  the 
radar  cross  section  of  three  dimensional  objects .  By  using 
the  natural  coordinate  system,  it  is  also  possible  to  analyze 
the  radar  cross  section  of  objects  of  composite  structure 
that  can  be  represented  by  nested  surfaces  of  arbitrary 
shapes . 

Extensive  benchmark  calculations  have  been  carried  out 
using  the  code.  Comparison  with  analytical  solutions  for 
scattering  off  simple  objects  showed  excellent  agreement  to 
less  than  1%  error  with  relatively  trivial  computational 
requirements  and  effort.  Results  of  calculated  radar  cross 
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section  for  ellipses,  rectangles  and  crosses  are  compared  to 
show  their  distinctly  different  signatures.  The  availability 
of  such  a  code  capable  of  efficient  and  accurate  calculation 
of  radar  cross  sections  can  be  expected  to  have  application 
to  defense  programs  in  three  important  areas.  The  first  is 
for  use  in  interactive  design  of  flight  vehicles  to  rapidly 
evaluate  different  design  options.  The  second  is  for  use  in 
interactive  design  of  detection  facilities  for  maximum 
sensitivity  and  discrimination  capabilities.  Thirdly,  with 
further  refinements  to  the  code  and  advances  in  computer 
hardware,  it  should  be  possible  to  implement  a  system  on 
manned  and  unmanned  flight  and  space  vehicles  with  capability 
for  on-line  processing  of  radar  signals  for  more  selective 
targeting  and  defensive  manuvers . 
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INTRODUCTION 


I . 

The  ability  to  calculate  accurately  the  radar  cross 
section  of  objects  of  arbitrary  shape  with  complex 
composite  structures  is  of  interest  to  defense  programs  for  a 
variety  of  reasons.  In  the  design  of  flight  vehicles  with 
minimum  radar  cross  sections,  it  is  desirable  to  be 
able  to  interactively  evaluate  the  design  options  by 
numerical  analysis  before  committing  to  hardware  modeling 
studies.  For  the  design  of  detection  systems,  it  is  equally 
important  to  be  able  to  determine  the  response  of  various 
types  of  threats  in  order  to  properly  configure  the 
detection  systems.  In  general,  both  monostatic  cross 
sections  of  backscattered  waves  and  multi-static  cross 
sections  as  seen  by  remote  receivers  are  of  interest. 

In  general,  the  radar  cross  section  of  an  object  can  be 
determined  from  the  solution  of  the  electromagnetic  wave 
equation  subject  to  the  boundary  conditions  at  the  surface  of 
the  scattering  object.  For  objects  of  very  simple  shape,  such 
as  spheres  and  infinite  cylinders,  analytical  solutions 
exist  from  which  the  differential  radar  cross  sections 
at  any  scattering  angle  can  be  determined.  In  the  more 
general  case,  the  solutions  can  be  found  only  by 
numerical  techniques.  Furthermore,  for  complex,  three 
dimensional  objects  with  composite  structures,  numerical 
solutions  by  a  brute  force  approach  can  be  taxing  for  even 
the  most  advanced  computers  and  they  are  not  well  suited 
for  interactive  or  field  usage.  Therefore,  there  is  the  need 
for  more  efficient  algorithms  for  the  accurate  solution  of 
Maxwell's  equations. 

In  recognition  of  such  needs,  we  have  developed  a 
significantly  more  efficient  approach  for  the  solution  of 
Maxwell's  equations  to  calculate  radar  cross  section  of 
three  dimensional  objects  of  arbitrary  shape  and  composition. 
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A  feature  of  the  approach  is  the  possibility  to  analyze  the 
response  from  nested  surfaces  of  composite  structures.  With 
support  from  the  Defense  Advanced  Research  Project  Agency 
(DARPA)  under  the  1988  DOD/SBIR  Phase  I  program, 
have  developed  a  formulation  of  this  approach  to 
demonstrate  its  capabilities. 

In  the  spirit  of  the  SBIR  program,  which  specifies  that 
the  objective  of  the  Phase  I  program  is  to  establish  insofar 
as  possible,  the  scientific  or  technical  merit  and 
feasibility  of  the  idea,  we  limited  our  objective  of  the 
six-month  Phase  I  effort  to  developing  the  algorithm  for 
calculating  the  cross  section  for  the  scattering  of 
electromagnetic  waves  from  infinitely  conducting  surfaces  of 
arbitrary  two-  dimensional  shapes.  The  formulation  can  then 
be  generalized  for  treating  scattering  from  fully  three- 
dimensional  objects  with  composite  structures  under  the 
more  extended  Phase  II  efforts.  The  results  of  the  Phase  I 
effort  is  fully  described  in  this  Phase  I  Final  Report. 
The  mathematical  concept  underlying  our  approach  is  described 
in  Section  II.  Details  of  the  formulation  and  the 
algorithms  for  solutions  can  be  found  in  the  Appendices. 
Examples  of  results  from  benchmark  calculations  are  presented 
in  Section  III.  Conclusions  regarding  the  demonstrated  and 
projected  capabilities  of  the  methodology  are  discussed  in 
Section  IV. 

II.  MATHEMATICAL  CONCEPT 

The  basic  problem  we  want  to  address  is  to  solve  the 
electromagnetic  wave  equation  in  infinite  space,  subject  to 
boundary  conditions  imposed  by  the  scattering  object  of 
arbitrary  shape  in  which  the  material  parameters  may  be 
inhomogeneous  and  anisotropic  with  sharp  discontinuities. 
It  is  intuitively  obvious  that  the  accuracy  of  the 
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solution  is  critically  dependant  on  the  proper 
representation  of  the  shape  of  the  scattering  object 
in  the  computational  grid  used  for  the  numerical 
solution.  For  objects  with  a  composite  structure,  such 
shapes  constitute  a  nested  sets  of  surfaces.  Therefore, 
accurate  representation  of  arbitrary  shapes  by  a  brute  force 
technique  would  require  the  use  of  a  very  fine  grid  which  can 
be  computationally  demanding. 


To  avoid  this  difficulty,  we  formulate  the 
electromagnetic  wave  scattering  problem  in  a  natural 
coordinate  system  that  conforms  to  the  arbitrary  shapes  of 
the  nested  series  of  surfaces.  Far  away  from  the  scattering 
object,  we  can  also  expect  the  effects  due  to  the  detailed 
features  of  the  scattering  surfaces  to  be  smoothed  out  and 
it  should  then  be  possible  to  solve  the  wave  equations  on  a 
more  conventional  grid.  To  minimize  the  computational 
efforts,  we  divide  the  infinite  space  into  two  regions 
separated  by  a  fictitious  boundary  that  fully  encloses  the 
scattering  object  (or  circle  for  a  two  dimensional  object) . 
Inside  the  fictitious  interface,  we  solve  the  wave  equations 
in  the  natural  coordinate  system.  Near  the  interface,  the 
natural  coordinate  system  gradually  conforms  to  the  shape  of 
the  interface.  Outside  of  the  interface,  the  wave  equation 
can  be  solved  in  a  conventional  coordinate  system  in  which 
the  functional  form  of  the  solution  is  known  and  the  boundary 
conditions  are  simple.  A  unique  solution  is  then  obtained  by 
matching  the  two  solutions  at  the  fictitious  interface.  By 
using  this  approach,  we  are  able  to  achieve  higher  accuracy 
with  less  computational  effort. 


Following  this  approach,  the  method  of  solution  consists 
of  several  steps.  The  first  step  is  to  generate  the  natural 
coordinate  system  in  the  multiply  connected  region  between 
the  scattering  object  and  the  fictitious  interface.  We  then 
rewrite  the  electromagnetic  wave  equation  in  this  new 
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coordinate  system.  Because  of  the  arbitrary  shape  of  the 
scattering  object,  the  boundary  conditions  are  complicated, 
which  will  be  governed  by  simple  boundary  conditions. 
However,  the  natural  coordinate  system  leads  to  simple 
boundary  conditions  and  the  wave  equation  can  be  easily 
solved  numerically  in  the  transformed  region.  Outside  the 
fictitious  interface,  the  solution  can  always  be  represented 
by  a  harmonic  series.  The  amplitudes  of  each  of  the  harmonics 
is  determined  by  matching  to  the  interior  solution  at  the 
fictitious  interface. 

For  the  Phase  I  effort,  we  restrict  ourselves  to  address 
the  two  dimensional  problem  of  calculating  the  scattering 
cross  section  for  infinite  rods  of  arbitrary  but  constant 
cross  sectional  shape.  We  chose  as  the  fictitious  interface 
an  infinite  cylinder  of  radius  a  that  fully  surrounds  the 
scattering  object  and  defines  a  doubly  connected  interior 
region  as  shown  in  Fig.  1.  We  shall  use  the  axis  of  the 
infinite  cylinder  as  the  origin  of  a  cylindrical  reference 
system  (r,(|»)  with  the  ignorable  coordinate  z  along  the  axis 

of  this  fictitious  cylinder. 


z 
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We  assume  there  is  a  plane  wave  of  specified  amplitude, 
polarization  and  angle  of  incidence  given  by  the  wave  vector 
ko  propagating  through  the  space  and  interacting  with  the 
scattering  object.  The  fields  outside  the  boundary  of  the 
cylinder  may  be  written  as  the  sum  of  the  incident  wave  in 
terms  of  the  components  of  the  wave  vector  ky,  kj  and  the 
wave  frequency  60;  and  the  scattered  wave,  which  we  shall 

express  as  a  superposition  of  cylindrical  waves  of  a  known 
analytical  form  obtained  from  the  solution  of  the  wave 
equation  incident  on  an  ideal  cylinder.  The  amplitude  of  this 
scattered  wave  is  unknown  and  it  contains  information  on  the 
geometry  and  material  characteristics  of  the  scattering 
object.  In  terms  of  the  (r,<())  coordinate  system,  the  incident 

field  is  represented  by 

£i^^-i(r.k«,-<ot)(Ej|e^+EiLe^)  (2.1) 

where  E^H  and  E^i  are  the  parallel  and  perpendicular 
components  of  the  incident  field,  and  ko=  (|X€)^^2  co/c  where  CO 
is  the  frequency  of  the  incident  wave. 

The  scattered  field  at  a  point  P  in  the  region  outside 
the  matching  interface,  can  be  written  as 


oo 

na-oe 


(2.2) 


where  (k^r)  is  the  cylindrical  Hankel  function  of  the  first 
kind  of  order  n  and  argxament  k^r.  The  constants  are  complex 
numbers  that  depend  on  the  geometry  and  the  material 
characteristics  of  the  scattering  object.  In  the  far-zone, 
which  is  when  the  conditions  kor  >>(koa)2  and  r  >>  X,©  are 

satisfied,  the  E-field  can  be  approximated  by  a  large 
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argument  expansion  for  the  Hankel  function  and  its  analytical 
form  can  be  written  as 


i  (k^r  -It  /4) 

— ==— T((|)) 

•/m 


(2.3) 


where 

T  (<)>)=  X  A„  e‘  ^ 


(2.4) 


The  scattering  cross  section  per  unit  length  can  be  obtained 
as 


CF(<j>)  =  ^  I  T(0)  P  (2.5) 

The  problem  is  now  reduced  to  solving  the  wave  equation 
inside  the  cylindrical  boundary  of  radius  r©,  and  matching 
the  solution  (with  a  scaling  factor)  to  the  sum  of  the 
incident  plane  wave,  (2.1),  and  the  scattered  cylindrical 
wave,  (2.2),  at  the  fictitious  cylindrical  surface  where 
there  is  no  source. 

In  the  more  general  case,  we  can  model  the  scattering 
object  as  a  set  of  nested  sharp  interfaces  across  which  the 
properties  of  the  media  changes  discontinuously .  A  natural 
set  of  coordinates  is  one  that  conforms  to  this  set  of  nested 
interfaces.  Each  surface  is  represented  by  setting  one  of  the 
coordinates  as  a  constant.  We  first  identify  the  boundaries 
of  the  sharp  interfaces  of  abruptly  changing  properties  in 
the  scattered  object  and  assume  that  the  material  properties 
between  the  sharp  interfaces  is  inhomogeneous  and  anisotropic 
but  varies  in  a  continuous  manner.  We  define  T|  as  a 

coordinate  function  which  is  constant  on  each  of  the 
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Fig.  2  Transformation  of  Physical  Plane  into  Nat' 
Coordinates . 


interfaces.  A  second  coordinate  function  ^ is  chosen  to 
intersect  the  resurfaces  on  the  interfaces,  and  4  =  4>  on  the 
cylindrical  boundary.  We  require  both  r|  and  4  satisfy 
Laplace's  equation 

V^T1  =  0  (2.6) 

and 

V^4  =  0  (2.7) 

subject  to  Dirichlet  boundary  conditions.  In  general,  the 
generated  coordinate  system  will  be  orthogonal  only  if  we 
require 

=  0  (2.8) 


dn 
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However,  this  requirement  would  not  allow  for  arbitrary 
spacing  of  the  natural  coordinate  lines  around  the  boundary. 
In  the  general  case  when  this  condition  is  not  satisfied,  the 
coordinate  system  will  not  be  orthogonal.  The  system  of 
governing  equations  must  then  be  transformed  directly 
through  implicit  partial  differentiation  rather  than  by  use 
of  scale  factors  and  differential  operators  as  is  possible 
with  orthogonal  curvilinear  systems.  However,  the  advantages 
of  allowing  for  arbitrary  spacing  far  outweighs  the 
complications  of  the  transformation  procedure  which  is 
carried  out  analytically  only  once.  Therefore,  we  shall  use 
the  more  general  non-orthogonal  curvilinear  coordinate 
system.  An  example  of  the  natural  coordinate  system  in  the 
physical  plane  f  an  irregular  object  is  shown  in  Fig.  2. 

The  effect  of  the  transformation  in  this  two  dimensional 
case  is  tc  conv'srt  a  doubly  connected  region  bounded  by  two 
closed  contours  of  arbitrary  shape  in  the  physical  plane  into 
a  rectangular  region  in  the  transformed  plane  in  which  two 
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of  the  opposite  boundaries  are  coincident  in  the  physical 
plane.  As  has  been  noted,  the  transformation  must  be 
carried  out  by  taking  the  partial  derivatives  according  to 
the  chain  rule.  In  general,  the  resulting  equation  will  be 
considerably  more  complicated  than  the  original  form. 
However,  the  boundary  conditions  are  now  specified  along 
straight  lines  and  the  equations  are  to  be  solved  in  a  simple 
rectangular  system. 

The  next  step  is  to  represent  the  Maxwell  equations  in 
the  natural  coordinate  system  and  determine  the  wave  equation 
for  either  the  electric  or  the  magnetic  fields  in  the 
interior  of  the  fictitious  cylinder.  We  will  also  base  our 
formulation  on  the  tensorial  representation  of  the  equations 
and  vector  fields,  since  all  the  information  about  the 
geometry  of  the  system  can  be  contained  in  the  components  of 
the  metric  tensor  which  simplifies  the  formulation  of  the 
problem  and  the  numerical  algorithms,  as  a  consequence. 

For  our  present  purpose,  we  shall  consider  the 
simplified  problem  of  a  two  dimensional  scattering  object 
that  is  a  perfect  conducting  infinite  cylinder  of  arbitrary 
cross  sectional  shape  and  solve  for  the  wave  equation  in  the 
multiply  connected  region  bounded  by  the  scattering  object 
and  the  minimum  cylinder  that  completely  encapsulate  this 
object . 

The  Maxwell  equations  in  the  source  free  space  are: 


V*E=0 

(2.9) 

_  1  ae 

V  X  E  - - ^ 

C  at 

(2.10) 

V  .  B  =0 

(2.11) 

_  _  1  aE 

V  X  B  =  -3- 

c  at 

(2.12) 
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where  K  ■  +  ■*  ,  and  E®  are  the  incident  and  scattered 
component  of  the  electric  field  respectively.  The  same 
notation  stands  also  for  the  magnetic  field  in  the  region. 


The  boundary  conditions  at  the  surface  of  the  perfect 
conducting  scatter  are  given  by 


n 

• 

E 

=  4ll  q. 

(2.13) 

n 

X 

E 

=  0 

(2.14) 

n 

• 

B 

=  0 

(2.15) 

B 

4n  , 

n 

X 

=  T 

(2.16) 

where  en  is  the  unit  vector  normal  to  the  scatter  surface, 
and  directed  outward,  is  the  surface  charge  density,  and  Js 
is  the  surface  current  density.  From  eqs.(2.10)  and  (2.12),  a 
wave  equation  for  the  electric  field  can  be  written  as 

VxVxE  =  -4^  (2.17) 

We  shall  represent  this  equation  in  terms  of  the  natural 
coordinates  by  making  use  of  the  conformal  transformation 
described  in  Appendix  B.  In  order  to  simplify  the  analysis, 
and  to  be  able  to  present  bench  mark  cases,  we  will  assume 
that  the  E-field  of  the  incident  plane  wave  is  linearly 
polarized  in  the  direction  parallel  to  an  infinitely  long 
circular  cylinder  which  is  a  perfect  conductor.  Therefore, 
the  scattered  electric  field  will  also  be  polarized  f^irallel 
to  the  cylinder  axis.  The  wave  equation  given  by  (2.17) 
becomes  a  scalar  equation  for  E®,  which  is  the  scattering 
field  in  the  physical  domain  (x,y),  and  it  can  be  written  as 

vX  +  i4e;;.=o  ,,,o. 


with  the  boundary  conditions 


[  Ejn  (*.y)  Ibody  —  [  E  (x,y)  Ibody  |2  i  gj 

t  (*.y)  IRo  =  [  Epu,  ((|),r  )  ]r^  (2  20) 

where  Eout  is  the  scattering  field  outside  the  matching 
boundary,  of  radius  Rq,  and  is  the  incident  plane  wave.  We 
also  require  that  the  normal  derivative  of  Ein  be  continuous 
at  Rq,  that  is. 


Ro 

dn 

Ro 


(2.21) 


When  transformed  to  the  natural  coordinate  system,  by  making 
use  of  Eq.B.13  in  appendix  B,  the  above  wave  equation  (2.18) 
becomes 


aE;;-2pE^^+YE^T,  +  J^ko  E  =  0 

(2.22) 

where  E  is  a  function  of  (^,  11  ),  and  E§^ ,  E§,,  and  £,,11  are  its 
partial  derivatives.  The  metric  elements  and  the  Jacobian,  J, 

the  transformation  are  calculated  in  Appendix  B  and  can 
be  written  as 

2  2 

(2.23) 

(2.24) 

-< 

II 

+ 

.A*fJ 

(2.25) 

(2.26) 
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The  transformed  boundary  conditions  are 


[  (i  n)  111,  =  [  E‘ (4, -n  )]ti,  (2.27) 

[  (;,T1)  ]tio  =  t  E^ut  ]r^  (2.28) 


'  aE:„,(<j>,r)' 

L  J 

Tio 

dr 

where  >  <t>  when  T|  — >  Tjo  .  The  coefficients  O,  P  ,  y  and  the 
Jacobian  are  calculated  during  the  generation  of  the  natural 
coordinate  system. 

After  the  mesh  is  generated,  the  Helmholtz  equation  in 
natural  coordinates,  (2.22),  is  to  be  solved  by  using  finite 
difference  techniqes.  Since  the  factors  (X,p,yand  J  are 
function  of  the  two  coordinates  ^  and  T),  the  Helmholtz 
equation  has  variable  coefficients.  However,  this  presents 
no  difficulty  in  the  numerical  solution. 

The  Helmholtz  equation  can  be  approximated  as  a  set  of 
linear  algebraic  equations  by  the  use  of  finite  differences. 
The  matrix  resulting  from  the  finite  difference  equations  is 
sparse  and  well  conditioned.  As  a  result,  standard  sparse 
matrix  subroutines  that  perform  Gaussian  elimination  can  be 
used  to  solve  the  set  of  linear  equations.  The  use  of  a 
sparse  matrix  solver  makes  the  overall  method  very  efficient 
with  respect  to  computational  resources.  Test  cases 
demonstrating  the  accuracy  of  the  solution  are  presented  in 
the  following  chapter. 
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III.  EXAMPLES  OF  BENCH  MARK  CALCULATIONS 

To  carry  out  the  cross  section  calculations  based  on  our 
approach,  we  had  to  develop  two  computational  modules.  The 
first  module  is  for  generation  of  the  natural  coordinate 
system  in  the  region  between  the  scattering  object  and  the 
fictitious  interface  and  transform  of  the  wave  equation  in 
this  region  into  the  natural  coordinate  system.  The 
fomnulation  for  the  determination  of  the  natural  coordinate 
system  is  described  in  detail  in  Appendix  A.  The  algorithm 
for  the  transformation  of  the  wave  equation  into  the  natural 
coordinate  system  is  described  in  Appendix  B.  The  numerical 
technique  for  the  solution  of  the  wave  equation  in  the 
natural  coordinate  system  and  matching  of  the  solution  at  the 
fictitious  interface  to  the  series  solution  outside  the 
interface  is  described  in  Appendix  C.  In  this  section,  we 
present  examples  of  calculated  results  and  benchmark  cases  to 
illustrate  the  application  of  the  code  to  calculate  the  radar 
cross  section  of  two  dimensional  objects  of  arbitrary  shape. 

A.  Generation  of  the  Natural  Coordinate  System 

The  objective  for  defining  a  natural  coordinate  system 
is  to  transform  the  governing  equations  from  doubly  connected 
region  R  bounded  by  contours  (surfaces  for  the  fully  three 
dimensional  case)  of  arbitrary  shape  into  a  rectangular 
region  R' .  The  effect  is  illustrated  in  Fig.  3  that  relates 
the  region  between  a  scattering  object  and  an  enclosing 
circle  in  the  (r,0)  or  (x,y)  physical  space  to  a  rectangle  in 
(T|,^ )  natural  coordinate  space.  The  natural  coordinates  (Tj.^) 
are  obtained  from  solutions  of  Laplace's  equations  for  T|  and  ^ 
with  Dirichlet  boundary  conditions:  T]  =  constant  -  T]i  on  Ci  , 
T|  *  constant  -  1^2  on  C2;  ^(x,y)  is  a  multiple  valued  solution 
with  a  branch  of  ^(x,y)  specified  (but  not  constant)  on  Ci  and 
C2  .  The  curve  Ci  on  the  physical  plane  transforms  to  the 
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Fig.  3  Mapping  of  the  coordinate  transformation 
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lower  boundary,  Ci' ,  of  the  transformed  plane.  Similarly,  C2 
transforms  to  C2 ' ,  etc.  The  right  and  left  boundaries  of  the 
rectangular  transformed  plane,  C3’  and  C4 '  ,  are  coincident  in 
the  physical  plane.  The  curve  which  transforms  to  these 
boundaries  connects  Ci  and  C2  and  determines  a  branch  cut  for 
the  multiple  valued  function  ^(x,y).  Thus  the  functions  and 
all  derivatives  are  continuous  across  this  cut. 

The  natural  coordinate  system  so  generated  has  a 
constant  T|-line  coincident  with  each  boundary  in  the  physical 
plane.  The  ^-line  may  be  spaced  in  any  manner  desired  around 
the  boundaries  by  specification  of  the  ^  boundary  conditions, 
or  equivalently  by  specification  of  (x,y)  at  the  equi-spaced 
^-points  on  the  tlx  and  T\2  lines  of  the  transformed  plane. 
Control  of  the  spacing  of  the  'H-lines  may  be  exercised  by 
varying  the  elliptic  system  of  which  T|  and  ^  are  solutions. 

To  test  the  module,  we  have  carried  out  the  calculation 
for  the  generation  of  the  natural  coordinate  grids  for  a 
number  of  two  dimensional  objects  with  geometric  cross 
sections  of  circles,  ellipses,  rectangles  and  crosses.  To 
speed  convergence,  we  have  found  founding  of  sharp  corners  as 
shown  to  be  helpful  but  not  necessary.  Exaunples  for  each  of 
the  four  cases  are  shown  in  Fig.  4,  5,  6  and  7,  respectively. 
All  computations  were  carried  out  on  an  IBM  3090.  Using  the 
algorithm,  convergence  to  less  than  1%  deviation  is 
achieved  in  only  a  few  iterations  for  the  elliptic  case  as 
shown  in  Fig.  8a  and  8b  where  we  have  plotted  the  normalized 
changes  in  the  x  and  y  coordinates  of  the  grid  for  the 
natural  coordinate  system  as  a  function  of  the  number  of 
i-erations.  Convergence  to  0.2%  is  realized  typically  in  30 
iterations  and  the  accuracy  continues  to  improve  with  further 
iterations.  Similar  behavior  are  found  for  the  other  cases 
where  we  have  carried  out  bench  mark  calculations  to 
determine  the  number  of  iterations  and  computer  processor 
unit  (CPU)  time  in  seconds  on  the  IBM  3090  as  a  function  of 
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convergence  tolerance  specified.  The  results  are  shown  in 
Fig.  9a,  9b  and  9c  for  the  case  of  the  ellipse,  cross  and 
rectangle,  respectively.  In  all  cases,  the  CPU  time  is  less 
than  a  few  seconds  for  convergence  tolerance  of  less  than  one 
part  per  million.  At  this  point  in  the  development  of  the 
code,  we  have  made  no  effort  to  optimize  the  code  and  even 
so,  the  computational  requirement  for  the  generation  of  the 
natural  coordinate  system  is  quite  insignificant. 

B.  Verification  of  the  Helmholtz  Solver 

Making  use  of  the  natural  coordinate  system,  the  next 
task  is  to  transform  the  governing  equations  into  the  natural 
coordinate  system  and  solve  the  system  of  equations 
numerically  subject  to  the  transformed  boundary  conditions. 
To  test  the  algorithm,  we  carried  out  a  benchmark  calculation 
of  a  cylindrical  stationary  wave  subject  to  simple  boundary 
conditions  for  which  an  exact  analytical  solution  is 
available  for  comparison. 

The  governing  equation  is  the  Helmholtz  equation  which 
has  the  form 


=  0 


«x 


In  cylindrical  coordinates  this  becomes 

idr^E„'i  1  a^E, 


— r 


^  a(j)' 


r  +  k2Eex  =  0 


r  9r 

Assuming  Eex  the  form 

E.x  =  E(r)cos<^ 

The  solution  of  the  equation  is 

E.x  =  {AJi(kr)  +  BYi(kr)}cos(j> 


(3.1) 


(3.2) 


(3.3) 

(3.4) 
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implies 


and 


_ -Yi(k) 

Ji(k)Yi(kro)-Ji(kro)Yi(k) 


(3.7) 


B  = 


Ji(k) 

Ji(k)Yi(kro)-Ji(kro)Yi(k) 


(3.8) 


which  gives  the  exact  solution. 


The  numerical  solution  is  obtained  by  transforming  the 
Helmholtz  equation  into  the  natural  coordinate  system.  The 
resulting  equation  in  the  (T|,^)  space  has  the  form 


aE^^-2PE^^  +  7E^^  +  k¥E  =  0 


(3.9) 


with  the  boundary  conditions 

1 11-1  =  0  and  E(^,ti)  1  =  cos  % 


(3.10) 


The  equations  are  solved  numerically  and  the  results  are 
compared  with  the  exact  analytical  solution  as  shown  in  Fig. 
10a,  10b  and  10c  in  which  we  compared  the  solutions  along  a 
radial  cut  at  45-degree  angle  for  three  cases  by  taking  20, 
40  and  60  radial  grids.  The  points  represented  by  the  black 
diamonds  are  the  numerical  results  while  the  open  squares 
represent  the  exact  analytical  values  calculated  at  the 
corresponding  grid  position.  The  improved  agreement  of  the 
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HELMHOLTZ  EQUATION  (  ROUT  s  10  ) 
(  PHI  =  45.0  DEGREES  :  JMAX  =  60  ) 
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Comparison  of  analytical  and  numerical  solutions  of 
Helmholtz  equation  at  a  constant  radius  for  (a)  20, 
(b)  90,  and  (c)  60  radical  grids 


numerical  solution  with  the  analytical  solution  with 
refinement  to  the  radial  grid  is  readily  apparent.  Excellent 
agreement  is  realized  by  taking  60  mesh  points.  The  agreement 
between  the  solutions  as  a  function  of  the  angle  at  a 
constant  radius  for  the  three  cases  are  shown  in  Figs.  11a, 
11b  and  11c.  It  can  be  seen  that  the  trend  is  similar. 

Based  on  these  results,  we  concluded  that  by  chosing  a 
fine  enough  mesh,  the  numerical  solver  will  match  to  the 
exact  solution.  This  is  due  to  the  fact  that  the  (x,y)  points 
generated  have  an  error  of  0(h2)  where  h  *  the  step  size  in 
the  radial  direction.  So  as  the  number  of  T|-lines  is 
increased  h  will  decrease  and  the  computed  solution  will 
become  a  better  approximation  to  the  exact  solution.  When  a 
mesh  of  ^max*  53  and  Tlniax  “  42  is  used  a  mean  difference 
between  the  exact  solution  and  the  computed  solution 
normalized  with  respect  to  the  exact  solution  of  4.5%  while 
with  a  mesh  of  ^max=  53  and  Tjinax  =  62  the  mean  difference 
between  the  exact  solution  and  the  computed  solution 
normalized  with  respect  to  the  exact  solution  was  found  to  be 
2.6%.  This  implies  that  as  the  mesh  is  refined  the  computed 
solution  is  converging  to  the  exact  solution. 

It  was  also  observed  that  the  n\jmber  of  ^-lines  does  not 
have  as  great  an  effect  on  convergence  as  the  number  of  T)- 
lines.  In  going  from  a  mesh  of  ^ax  =  97  and  Tlmax  =  21  to  a 
mesh  of  ^ax  *  53  and  Tlmax  “  42  the  doubling  of  the  number  of 
T|-lines  overcomes  any  decrease  in  convergence  from  decreasing 
the  number  of  ^-lines. 

C.  Results  of  Cross  Section  Calculations 

To  establish  the  validity  of  the  code,  we  compare  the 
results  of  the  numerical  calculations  to  the  analytical 
results  for  scattering  from  a  perfect  conducting  infinite 
cylinder  of  circular  geometric  cross  section.  An  analytical 
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Helmholtz  equation  at  a  constant  radius  for  (a)  20, 
(b)  90,  and  (c)  60  radical  grids 


solution  for  the  case  when  the  incident  electric  field  is 
linearly  polarized  (parallel  to  the  axis  of  the  cylinder)  is 
given  in  Appendix  D,  which  is  of  the  form 


E*  =  5^  (-i)"  £„  A„  (koT)  cos  (n<|)) 

n=4) 


with 


„  /  1  f or  n  =  0 

£n  “1  2  for 1,2,3,... 


(3.11) 


where  An  is  the  coefficient  of  the  n'^^  harmonic  in  the  series 
representation  of  the  analytical  solution.  For  waves 
scattered  by  an  infinite  cylinder,  the  magnitude  of  An  as  a 
function  of  n  for  the  first  five  harmonics  are  shown  in  Fig. 
12a.  It  can  be  seen  that  for  low  wave  frequencies  (or  ka 
values)  the  magnitude  of  An  falls  off  very  rapidly  with 
increasing  n  such  that  the  contributions  from  high  order 
modes  will  be  negligible.  Fig.  12b  shows  that  for  higher 
frequencies  the  contributions  from  high  order  modes  is  very 
significant.  In  matching  the  interior  solution  to  the 
external  analytical  solution  at  the  boundary,  we  need  to 
evaluate  many  high  order  terms  for  higher  frequency  incident 
waves . 


The  total  electric  field  as  a  function  of  the  angle 
around  the  scattering  object  for  waves  incident  along  <j>o  =  0  is 
shown  plotted  in  Fig.  13.  The  corresponding  differential 
scattering  cross  section  as  a  function  of  the  scattered  angle 
is  shown  in  Fig.  14  for  the  case  k  =»  1 .  In  general,  for  high 
wave  numbers,  the  cross  section  is  maximum  in  the  forward 
direction.  The  forward  ((|)  =  n)  and  backward  (<|>  *  0)  scattering 
cross  section  as  a  function  of  the  wave  number  is  shown  in 
Figs.  15  and  16,  respectively.  These  results  calculated  from 
the  analytical  expression  provide  a  reference  for  evaluating 
the  accuracy  of  our  ntimerical  solution. 
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The  differential  scattering  cross  section  for  the 
infinite  cylinder  calculated  with  the  code  was  used  as  a 
bench  mark.  The  results  for  the  case  of  an  incident  wave 
with  wave  number  k  =  1  is  shown  in  Fig.  17  in  which  the 
differential  scattering  cross  section  as  a  function  of  the 
scattered  angle  is  given  by  the  curve  on  top  and  the  natural 
coordinate  grid  used  for  the  calculation  is  shown  at  the 
bottom.  The  results  calculated  with  the  code  is  compared  with 
the  cross  section  obtained  from  the  analytical  solution  as 
shown  in  Fig.  18.  The  absolute  value  of  the  normalized  cross 
section  are  compared  in  the  top  figure  and  their  differences 
as  a  function  of  the  scattered  angle  are  plotted  in  the 
bottom  figure.  It  can  be  seen  that  the  agreement  is  very  good 
and  the  deviation  is  less  than  1%  at  all  scattered  angles. 

To  further  verify  the  accuracy  of  the  numerical 
procedure,  we  evaluated  the  coefficient  of  the  individual 
harmonics  at  the  matching  interface  for  comparison  with  the 
coefficients  determined  from  the  analytical  solution.  The 
results  are  shown  in  Fig.  19  for  the  case  ka  =  1 .  It  is  seen 
that  the  agreement  is  very  good  for  the  first  five  harmonics. 
Since  the  contributions  of  the  high  order  harmonics  is 
totally  negligible  in  this  case,  it  is  reasonable  to  expect 
the  high  degree  of  agreement  between  the  two  results. 

To  exercise  the  code,  we  calculated  the  differential 
scattering  cross  sections  for  a  variety  of  shapes  by 
systematically  varying  the  wave  number,  the  angle  of 
incidence  and  the  number  of  modes  included.  In  general,  for 
ka  <  5,  we  found  that  higher  order  terms  has  negligible 
effect  on  the  solution  as  was  to  be  expected.  The  cases 
calculated  are  summarized  in  Table  1.  The  results  are  shown 
in  Figs.  20  through  30  and  a  variety  of  trends  can  be 
observed,  some  are  intuitively  obvious  and  many  are  quite 
interesting . 
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Fig.  19  Comparison  of  the  harmonic  coefficient 

for  the  analytical  and  numerical  solutions 
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Table  1.  Summary  of  Cross  Section  Calculation. 


Wave  No.  of  Angle  of 

Fig.  No.  Shape  Number  Modes  Incidence 

20  Ellipse  0.5  6  0* 


21  Ellipse  1.0  6  0* 


22  Ellipse  3.0  6  0* 


23  Ellipse  5.0  6  0* 


24  Short  1.0  4  O' 

Rectangle 

25  Long  1.0  4  0' 

Rectangle 

26  Short  Arm  1.0  4  0' 

Cross 

27  Short  Arm  1.0  4  45* 

Cross 

28  Short  Arm  1.0  4  90* 

Cross 

2  9  Short  Arm  1.0  4  Varied 

Cross 

30  Long  Arm  1.0  4  0’ 
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Fig.  21  Cross  section  of  a  1x2  elliptic  cylinder  for  )<a=l  ' 
and  n  =5  modes  . 
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Fig.  22  Cross  section  of  a  1x2  elliptic  cylinder  for  ka-3 
and  n=6  modes  . 


Data  from  ‘'ELI  A05.XSEC 
WAVN  =  5  .  NCOEF  =6 


Data  from  ''HC"!  A01.XSEC 


Fig . 


24  Cross  section  of  a  1x2  rectangular  cylinder  for 
ka=l  and  n=4. 


Fig.  29  Cross  section  ac  (Me  of  the  geometrical  cross  as  a 
function  of  incident  angle  d> 
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Comparing  the  differential  scattering  cross  section  of 
an  ellipse  shown  in  Fig.  21  with  that  of  the  circle  shown 
previously  in  Fig.  17  for  the  same  wave  number,  we  note  the 
appearance  of  side  lobes  at  right  angle  associated  with  the 
eccentricity  of  the  shape.  The  dependence  on  the  wave 
frequency  for  scattering  by  the  same  elliptical  object  under 
same  angle  of  incidence  can  be  seen  by  comparing  Figs.  20, 
21,  22  and  23.  When  the  wavelength  of  the  incident  wave  is 
long  compared  to  the  scale  length  of  the  scattering  object, 
the  differential  scattering  cross  section  of  the  ellipse 
resembles  that  of  a  circle  for  scattering  of  higher  frequency 
of  waves.  As  the  wavelength  becomes  comparable  with  the  scale 
length,  the  forward  scattering  cross  section  rises 
dramatically  by  an  order  of  magnitude  and  side  lobes  appear 
at  right  angle. 

The  effect  of  the  length  of  the  scattering  object  along 
the  direction  of  incidence  can  be  seen  by  comparing  Figs  24 
and  25  for  two  rectangles  of  different  length.  The  magnitude 
of  the  forward  scattering  cross  section  scales  approximately 
linearly  with  the  length  of  the  object  while  the  percent 
increase  in  the  backward  scattering  cross  section  is  much 
more  pronounced.  Referring  back  to  the  cross  section  of  the 
ellipse  shown  in  Fig.  21  which  has  similar  aspect  ratio  as 
the  rectangle  in  Fig.  24,  we  note  that  in  fact,  the 
characteristics  of  their  differential  scattering  cross 
sections  are  quite  different. 

The  relative  effects  of  pretruding  arms  on  the 
scattering  cross  section  compared  to  length  can  be  deduced  by 
comparing  Fig.  25  for  the  long  rectangle  with  Fig.  26  for  a 
short  arm  cross  and  Fig.  30  for  a  long  arm  cross.  The  most 
notable  feature  is  that  increasing  the  length  of  the  body 
results  in  a  significantly  larger  increase  in  the  forward 
scattering  cross  section  than  the  effect  of  stubby  pretruding 
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side  arms.  Furthermore,  increasing  the  length  of  the  side 
arms  increases  the  forward  scattering  cross  section  but 
reduces  the  backward  scattering  cross  section,  which  is 
somewhat  counter  intuitive.  The  effect  of  the  angle  of 
incidence  on  the  scattering  cross  section  is  evaluated  by 
comparing  Figs.  26,  27,  and  28  as  the  angle  of  incidence  is 
increased  from  0  to  90  degrees  at  45  degree  intervals.  The 
general  trend  is  that  magnitude  of  the  cross  section 
increases  as  the  angle  of  incidence  is  increased  from  0  to  90 
degrees.  The  forward  scattering  cross  section  as  a  function 
of  angle  of  incidence  is  shown  plotted  in  Fig.  29.  It  should 
be  emphasized  that  this  figure  should  be  interpreted  with 
some  care  because  when  the  angle  of  incidence  in  180  degrees, 
the  "forward"  direction  is  actually  the  backward  direction. 

IV.  CONCLUSIONS 

The  present  work  was  carried  out  under  the  Small 
Business  Innovation  Research  Phase  I  program  which  is  limited 
to  six  months  in  duration  at  approximately  one-third  man-year 
level  of  effort.  Under  this  effort,  we  have  successfully 
demonstrated  the  concept  of  using  a  natural  coordinate  system 
that  conforms  to  objects  of  arbitrary  shape  for  the 
calculation  of  radar  cross  sections  of  such  objects.  The 
resulting  algorithm  is  shown  to  be  highly  versatile, 
efficient  and  accurate.  The  present  code  is  limited  to 
treating  two  dimensional  objects  of  arbitrary  shape  to 
establish  the  scientific  and  technical  merit  of  the  concept. 
The  formulation,  however,  is  completely  general  and  can  be 
extended  directly  to  the  full  three  dimensional  case. 

Making  use  of  the  code,  we  have  evaluated  the  radar 
cross  section  for  a  variety  of  shapes .  Comparison  with  known 
analytical  results  for  an  infinite  cylinder  showed  excellent 
agreement.  Reason  for  good  agreement  is  also  established.  To 
fully  test  the  code,  the  differentially  radar  cross  suction 
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for  a  variety  of  shapes  were  calculated  by  systematic 
variation  of  the  key  parameters.  Upon  examining  the  results, 
signatures  associated  with  the  different  shapes  can  readily 
be  identified,  suggesting  techniques  for  reducing  radar  cross 
sections . 

The  generalization  of  the  code  to  three  dimensions  is 
the  logical  next  step  that  can  lead  to  the  availability  of  a 
robust  and  efficient  code  for  the  calculation  of  radar  cross 
sections . 


5] 


APPENDIX 


APPENDIX  A 


GENERATION  OF  NATURAL  COORDINATE  SYSTEM 

A.l  Mathematical  Formulation. 

A  general  method  of  generating  natural  coordinate  system 
is  to  let  the  natural  coordinates  be  solutions  of  an  elliptic 
partial  differential  system  in  the  physical  plane,  with 
Dirichlet  boundary  conditions  on  all  boundaries.  The 
procedure  is  not  restricted  to  two  dimensions,  and  allows  the 
coordinate  tangential  to  the  boundary  to  be  distributed  quite 
easily  as  desired  along  the  boundary,  and  is  applicable  to 
all  mult iconnected  regions .  The  coordinate  system  so 
generated  is  not  necessarily  orthogonal,  but  orthogonality  is 
not  required,  and  its  lack  only  requires  that  the  governing 
equations  of  the  physical  system  must  be  transformed  directly 
through  implicit  partial  differentiation  rather  than  by  use 
of  the  scale  factors  and  differential  operators  developed  for 
orthogonal  curvilinear  systems.  An  orthogonal  system  cannot 
be  achieved  with  arbitrary  spacing  of  the  natural  coordinate 
lines  around  the  boundary,  and  the  capability  for  such 
arbitrry  spacing  is  of  more  importance  than  orthogonality. 

The  objective  is  to  transform  the  two  dimensional, 
double  connected  region,  R,  bounded  by  two  closed  countours 
of  arbitrary  shape  into  a  rectangular  region,  R'  ,  as  shown 
in  Fig.l.  The  general  transformation  from  the  physical  plane 
(x,y)  to  the  transformed  plane  (^TJ)  is  given  by 


(a.l) 
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Similarly,  the  inverse  transformation  of  (a.l)  is 


'xl  _rx(^,Tl)‘ 

.yJ  Ly(4.'n)- 


with  the  Jacobian  matrix 


Lny. 


which  yields  the  Jacobian  of  the  transformation 

J  =  x^y^  -  x^y§ 


(a. 2) 


(a. 3) 


(a. 4) 


The  partial  derivatives  of  a  generic  function  f  (x,y)  are 
transformed  by 


f  ^^.y)  /  ^x,y)  yti^4  ~ 

,  9(x,f)  a(x,y)  -  x^f,, 

'y "  a(iTi)  d(%.r\)  ■  J 


(a. 5) 


(a.  6) 


Since  the  basic  idea  of  the  transformation  is  to 
generate  transformation  functions  such  that  all  boundaries 
are  coincident  with  coordinate  lines,  the  natural  coordinates 
(4,^1)  are  taken  as  solutions  of  some  suitable  elliptic 
boundary  value  problem  with  one  of  these  coordinates  constant 
on  the  boundaries.  Using  Laplace's  equation  as  the  generating 
elliptic  system,  we  have 


^xx  +  ^yy  =  0 

(a. 7) 

Tlxx  +  nyy  =  0 

(a. 8) 
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with  Dirichlet  boundary  conditions,  T)  =  constant  =  Tj!  on  Ci  , 
11=  constant  =  TI2  on  C3  ;  ^  (x,y)  a  multiple  valued  solution 
with  a  branch  of  ^  (x,y)  specified  (but  not  constant)  on  Ci 
and  C3  .  The  curve  Ci  on  the  physical  plane  transforms  to  the 
upper  boundary,  ' ,  of  the  transformed  plane.  Similarly,  C3 
transforms  to  C3 ' ,  etc.  The  right  and  left  boundaries  of  the 
rectangular  transformed  plane,  C2'  and  C4 '  ,  are  coincident  in 
the  physical  plane.  The  curve  which  transforms  to  these 
boundaries  connects  Ci  and  C3  and  determines  a  branch  cut  for 
the  multiple  valued  function  4  (x,y).  Thus  the  functions  and 
all  derivatives  are  continous  across  this  cut. 

Now  since  we  wish  to  do  all  numerical  computation  in  the 
rectangular  transformed  plane,  it  is  necessary  to  interchange 
the  dependent  and  independent  variables  in  (a. 7)  and  (a, 8). 
Thus,  using  the  partial  derivative  transformations  and  the 
chain  rule,  equation  (b.l3)  in  Appendix  B  yields 

=  0  (a -9) 
-2py^^+ yy  1,^  =  0  (a.  10 


where  the  coefficients  a,  ^  and  y  correspond  to  the 
components  of  the  metric  tensor  of  the  coordinate 
transformation,  and  are  given  by 


a=x^  +  y5 

P  =  +  y^y^ 

y=x|  +  >i 


(a. 12) 
(a. 13) 
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with  the  transformed  boundary  conditions  given  by 


'xl  ni(^.Tli)' 

.yJ 


'xi  rgiC^-na)' 

.yJ  [92(4.^12). 


»  (4' '•12)  ®  ^3 


1 


(a. 14) 


This  system  is  a  quasi-linear  elliptic  system  with  Dirichlet 
boundary  conditions  for  the  physical  coordinates  in  the 
transformed  plane.  The  differential  equations  of  the  system 
(a. 9)  and  (a. 10)  are  considerably  more  complicated  than  those 
of  Eqs .  (a.  7)  and  (a.  8)  .  However,  in  the  first  case  the 
transformed  boundary  conditions  are  specified  on  straight 
boundaries,  and  the  computation  region  is  rectangular.  We 
have  thus  exchanged  a  problem  having  simple  equations  but 
complex  boundary  conditions  for  a  problem  having  complex 
equations  and  simple  boundary  conditions.  This  statement  also 
holds  for  all  partial  differential  equations  which  govern  the 
physical  system  solved  on  the  natural  coordinates . 

The  natural  coordinate  system  so  generated  has  a 
constant  Tl-line  coincident  with  each  boundary  in  the  physical 
plane.  The  ^-lines  may  be  spaced  in  any  manner  desired  around 
the  boundaries  by  specification  of  the  ^  boundary  conditions, 
or  eqpiivalent  by  specification  of  (x,y)  at  the  equi-spaced 
points  on  the  T|i  and  T|2  lines  of  the  transformed  plane. 
Control  of  the  spacing  of  the  T|-lines  may  be  exercised  by 
varying  the  elliptic  system  of  which  ^  and  T|  are  solutions. 

As  noted  above,  orthogonality  is  not  required  of  the 
coordinate  system  for  solution  of  a  system  of  partial 
differential  equations.  Indeed,  normal  derivatives  on  the 
boundaries  may  be  easily  represented  in  a  non-orthogonal 
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system  as  follows:  Let  y=  g(x),  i.e.,  y-g(x)-  f(x,y)*0,  be 
the  equation  of  some  boundary  in  cartesian  coordinates.  Then 


—  =  n.V  f  = 
9n 


- \  2  1/2  [  U  (9  Vn  +  ^  (g  i 

Jd-Kgry"^ 


(a. 15) 


where  g'  =  dy/dx.  All  derivatives  in  the  last  expression  can 
be  calculated  along  coordinate  lines  in  the  natural 
coordinate  system.  Thus  it  is  not  necessary  to  require 
orthogonality  at  the  body  surface  in  order  to  get  accurate 
representation  of  the  normal  derivative. 


A. 2  Numerical  Method. 


The  numerical  solution  of  the  system  of  Eqs.  (a. 9)  and 
(a. 10),  with  the  boundary  conditions  given  by  Eqs.  (a. 14),  is 
accomplished  by  the  use  of  a  nine-point  finite  difference 
approximation.  The  diagram  in  Fig.(a.l)  shows  the  nodal 
arrangement  for  the  approximation,  and  the  parameter  D  is 
the  separation  between  adjacent  nodes.  Now  using  central 
differences  to  aproximate  the  partial  derivatives  of  the 
variables  x  and  y  with  respect  to  ^  and  Tl  we  obtain  the 

following  expressions : 


^11 


Xg  —  2Xc  +  Xy^ 


XnE  ~  Xnw  ~  XsE  ^sw 

4a2 

xn-2xc  +  Xs 


Xg-Xw 

^  2A 

Xn-Xs 
■c.  * 

^  2A 


(a. 16) 
(a. 17) 
(a. 18) 
(a. 19) 
(a. 20) 
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And  substituting  the  above  expressions  into  Eqs.(a>9)  through 
(a. 14) ,  we  obtain 


a(XE  +  Xy^)  -  0.5p(X|^g  -  Xjgyy  -  Xgg  +  X3yy)  +  ^(Xi^  +  X3) 

2(a  +  y) 


1 


(a. 21) 


«(yE+yw)-o-5p(yNE-yNw-ysE-ysw)+7(yN+ys) 

2(a  +  y) 


(a. 22) 


where 

a  =  (Xn  -  Xq)^  +  (y^  -  Ysf 


(a. 23) 


p  =  (xe  -  xw)(xn  -  xs)  +  (yE  -  y  w)(yN  -  ys)  (a .  2  4 ) 

Y=(XE-Xw)^  +  (yE-yw)^  (a.  25) 

The  above  system  will  now  be  solved  on  the  (^>T|) 
rectangular  plane  with  =  if  for  i  =*  1,  2,  ...  ,M-1  and  T|j  - 

j,  for  j  =  1,  2,  ...  ,N.  There  are  4(M-1)  known  values 

corresponding  to  the  x  and  y  points  at  the  boundaries,  and 
therefore  there  are  2 (N-2) (M-1)  unknown  values  to  be 
determined.  We  will  use  Gauss-Seidel  iteration  scheme,  which 
is  a  direct  substitution  method  using  updated  values  as  they 
are  generated  in  an  iteration. 

The  values  of  x  and  y  on  the  left  and  right  boundaries 
of  the  transformed  plane,  C2  ’  and  C4 '  ,  are  unknown  but 
continuous  in  T|.  This  implies  that  the  values  of  xw,xse  and 
xne  at  i-M-1  are  also  the  values  i  -1.  These  relations  also 
hold  for  the  y-values  at  those  boundaries.  A  good  initial 
guess  is  needed  to  obtain  convergence  in  Gauss-Seidel  method. 
Since  the  x's  and  y's  in  the  radial  direction  will  fall 
between  the  values  at  the  boundaries  it  seems  that  a  weighted 
average  of  these  values  constitutes  an  appropriate  guess. 
Therefore,  the  values  for  Xij  and  yij  can  be  evaluated  from 
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(a. 26) 


^ij-  21 


N-i 

l.N-1 


ri1 


t^N-l  ) 


l^iN 


(a. 27) 


where  j=  2,3,...,N-1  and  i*l, 2, . . . ,M-1 . 


Comparison  of  the  difference  of  the  absolute  norm  of  the 
kth  and  k  -  1st  iterations  was  used  for  checking  convergence. 
The  II  X  II  is  calculated  by 


N-1M-1 

III  I 


IUII«=-i 


2U1 _ 

(M-1)(N-2) 


(a. 28) 


where  the  difference  in  the  norms  is  calculated  as 


DIff.  ||x||  =  |  ||xP>-||x||<''-^>|  (3,2! 

And  similarly  for  the  difference  in  the  y-norm.  If  the 
difference  in  the  norms  is  less  than  a  given  tolerance  the 
iteration  process  converges  and  the  final  values  for  x  and  y 
will  define  the  grid  points  corresponding  to  the  (^.Tl) 

coordinates . 

The  set  of  finite  difference  Eqs.(a.21)  and  (a. 22)  has 
been  successfully  solved  for  a  number  of  two  dimensional 
objects  with  various  geometric  cross  sections  corresponding 
to  circles,  ellipses,  rectangles  and  even  more  arbitrary 
geometries.  For  all  these  cases  the  outer  boundary  is  defined 
by  a  perfect  circle,  as  desired  by  the  matching  interface  in 
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the  scattering  problem.  The  plots  shown  are  of  constant 
lines  and  T|  =  constant  lines  in  the  physical  plane. 

The  specification  of  the  convergence  tolerance  was  found 
to  be  dependent  on  the  number  of  iterations.  The  number  of 
iterations  as  a  function  of  convergence  tolerance  has  the 
form  of  ce~^/^  (  where  x  =  convergence  tolerance  )  .  A  false 
convergence  was  found  when  the  niomber  of  iterations  was  below 
150.  This  implies  that  the  convergence  tolerance  must  be  at 
most  0(10“^).  Also  the  N-1  contour  is  required  to  be 
cylindrical  which  implies  that  the  deviation  in  the  radius 
must  be  small.  A  convergence  tolerance  of  0(10“^)  satisfied 
both  of  the  above  criterion. 

It  was  observed  that  the  number  of  points  per  contour 
had  a  more  significant  effect  on  the  cost  of  generation  than 
did  the  number  of  contours.  The  doubling  of  the  number  of 
contours  increased  the  number  of  iterations  by  less  than  5%. 
The  doubling  of  the  number  of  points  per  contour  doubled  the 
number  of  iterations  recjuired.  It  was  also  observed  that  the 
greater  the  number  of. contours  the  better  the  resolution  near 
the  inner  boundary. 
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APPENDIX  B 


Transformation  Algorithm 

I.  Derivatives  and  Vectors  in  the  Transformed  Plane 

Given  that  x=x  (^.^ )  and  y=y(4f^)  we  obtain  a  set  of 
relations  in  the  transformed  (x,h)  plane. Let  f  =  f(x,y)  -  a 
twice  differentiable  scalar  function  of  x  and  y. 


(  2  2  ^ 

x|  2*5/5  n 

V  ^  yri  / 
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Hi 

B=  y?, 

yiTT,/ 

detM=detA-  |j|  =  IjI'* 

By  inverting  the  5x5  matrix  M  we  obtain  the  relations  between 
the  derivatives  of  f  in  the  x-y  plane  with  respect  to  the 
derivatives  in  the  x-h  plane. 


(b.2) 


where 

IjI'C  d 

lo  0  0  Ijpgy 

^  Vn  2y^yr,  yf 
c=  -Xt,yn  -x^y^ 

,  x^  -2x5X^  x| 

D  =  rx^a-y^b  -x^a+y^b" 
-^o+y^6  x^c-y^d 

^x^o-y^g  -x^o+y^g> 

where 

a=yfy,,-2y5y,y5,  +  y5y« 
‘>=y|>s,-2y5y,*?.i  +  y^« 
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0  =  -(X|y,  +  x,y5)y5„  +  x^y^y,,, 

cl  =  x^yi,xu  -(x?y,  +  x,y|)x5,,  +  x^y^x^, 
9  =  x^y«  -  2x5x^y4,  +  xfy^^ 

9  “  Xt)X^^  —  2x^x^x^^  +  x^x,^^ 


From  (1)  and  (2)  we  obtain  the  derivative  transformations 
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The  Laplacian  is 


V^  = 


[(ox^^  -  2px^^  +  Yx^^)(y^f^  -  y„f^)] 


(gy^^  -  ->-  Yy^inXy^  " 

IJP 


p2.  +  (y^n+^^ 

I  ||2 


where 


2  2 
a  =  x;  +  y; 

P  =  +  y^yt, 

Y  =  xf  +  yf 

Dx  =  ox^^  -  2Px^^  +  yx^^ 

Dy  =  ay^^-2py^^+Yy^^ 
yeDx  -  x^Dy 


x^Dy-y^Dx 


In  particular  the  derivatives  of  4(x,y)  and  Ti(x,y)  are 


4xx  = 


ax^-byri 


(b.8) 


(b.9) 


(b.10) 


(b.11) 
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-  gy^ 

IjP 

2.  -  (t>+g)yii  X 

Ijp  "IJP 

Similarly 


(b.12) 


2  2 

If  V  4  =  0  V  T]  =  0  -hen  o  =  0  and  T  =  0  and  (9)  is  transformed 
into 


,2. 

IJP 


(b.13) 


II.  Unit  Tangent  and  Unit  Normal  Vectors  in  the  (x,h)  plane 


utilizing 

(3) 

and  (4)  we  have 

Yr\  -y^i+x^i 

(b.14) 

n  — 

1 

which  is 

the 

unit 

vector  normal 

to  a  line  of  constant  h. 

Similarly 

the 

unit 

vector  normal 

to  a  line  of  constant  x  is 

given  by 

n<^  = 

V4 

1  ^  1  “■  1— 

(b.15) 

“  lY4l  Va 


The  unit  tangent  vectors  are  given  by 

<y  (b.16) 

Va 
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APPENDIX  C 


METHOD  OF  SOLUTION  FOR  FINDING  FOURIER  COEFFICIENTS 


Because  of  the  linearity  of  the  governing  equations,  we  can 
write  the  total  field,  Etotf  as  the  sum  of  an  incident  field  and 
a  scattered  field. 

^tot  ~  ^in 

p  _  -ikrco8{Oo) 


+<>o 

E.=  E 

m=-<« 


The  constants,  A^,  are  generally  complex.  are 
cylindrical  Hankel  Functions  of  the  first  kind  of  order  m: 

and  e'"’^  ■  cos(m<)))  +  i  sin(m<j>) 

On  the  inside  boundary 

Etot  *  0^  therefore,  E,  =  -Ein 

The  scattered  field  can  be  broken  down  into  a  real  and 
imaginary  part : 


(E,)r  =  X  [a„,J„-b^YJcos(m(|))-[a„Y,„  +  bmJjsin(rTV|>) 
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+0O 

(Es)i  =  X  tamYm  +  b„JJcos(m(t>)  +  [a„J„,-b^YJsin(m<j>) 

m  =  -oa 

It  is  desirable  to  combine  terms  so  that  the  sxammation  is 
over  positive  integers. 

(E3)R  =  [aoJo-boYJ 

+  X  -  bmYm  "  b^(-1)"’Y^]  COS(m<|)) 

m=1 

+00 

-  X  +  b„j„  -  b^C-I^J^l  sitKm*) 

m=1 

(E3),  =  [aoYo  +  bM 

+00 

+  X  +  a_„(-1)"’Y„  +  b„Jm  +  cosCmij)) 

m*1 

+00 

+  X  “  a_n,(-1)'”j^  -  b^Y^  +  b_n,(-1)"’Y^]  Sin(m(j>) 

m=1 

For  the  ideal  case  of  an  obstacle  having  a  circular  shape,  these 
equations  can  be  matched  -  equated  to  Fourier  expansions  of  the 
boundary  conditions. 

For  example,  at  r*l, 

(Ea)R  =  -cos[k  cos(<j)  -  ^o)J  s 
(Eg)!  =  +sin[k  cos(<j>  -  <|>o)]  *  g(<(>) 
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Both  f(<{>)  and  g(<|))  can  be  expanded  as  Fourier  series: 


■foo 

f(<j>)  =  Ao  +  ^  cos(m<|))  +  Bn,  sin(m<|)) 

m=1 

+00 

g((|>)  =  Co  +  ^  Cm  cos(m4>)  +  sin(m(j>) 

m=1 


Matching  mode  by  mode  gives: 


[Aol 

1 

C- 

o 

1 

-< 

o 

_ 1 

'ao' 

1 - 

o 

o 

1 _ 

- 1 

o 

_ 1 

>o. 

and  for  m=l  to  +oo 


■■kj„  -Y„  -(-i)X' 

’am' 

Bm 

-Y„  -K-irY^  -K-irJ™ 

a-m 

E 

O 

+Y„  -K-I^Y^  +(-1)"'j„ 

m 

1 

E 

O 

- » 

.+Jm  -Ym  -K-DX. 

L  HTIJ 

If  we  were  solving  the  problem  in  cylindrical  coordinates,  the 
unknown  constants  {am,  a_m,  bm,  b_n,]  could  be  easily  determined 
by  solving  the  simple  4x4  system  of  linear  equations. 

For  the  general  boundary  shape  such  that 

robs  *  ri(^) 


The  problem  is  more  complicated, 
on  rin. 


(E»)r  =  -co8[k  r|(^)  cos(<|>  -  <^o)l 
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(E,)i  =  +sin[k  ri((j))  cos((t>  -  <|)o)] 


However,  in  the  natural  coordinate  system,  (%,X\) ,  the  boundary 
ri  ((|))  corresponds  to  Tj*'! .  We  are  at  liberty  to  define  the 
following  expansion: 


+»» 

f(V)  =  ^  +  Am  cx)s(m4‘)  +  Bm  sin(m4‘) 

m*1 


+» 

g(^*)  =  Co  +  Cm  cos(m4*)  +  Dm  sin(m^*) 

m=1 


where 


2n  (^-1) 


The  constants  S’^,,  D^,  are  line/contour  integrals  over  the 

boundary  shape . 

®-g-  ^  =  ^  f  cos(kr(4*)  cos(<|>(0  - <l>o)) 

^  Jo 


since  r(^*)  and  <|>(V)  points  on  the  surface  in  the  pareuneter 

V- 


For  this  to  work,  we  need  to  find  functions  as  the  following 
form: 


+0O 

[E,(^,ii)]r  =  {am  cx)s(m^*)  -  sin(mV)] 

-bm  [Bm(4,Tl)  cos(m^*)  -  Bm(4.Ti)  sin(m4*)]} 
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4«o 

[E3(4.ti)]i  =  X  ^ cos(m|‘)  +  /C(^.T1)  sin(m^‘)] 

m»-M 

[bJ,(4.ti)  cos(m^‘)  -  Bm(iTi)  sin(m4 )]} 

In  which  the  general  functions 

replace  the  Bessel  functions  used  in  cylindrical 

coordinates.  We  can  pick  the  natural  coordinate  system  in  such  a 
way  that,  at  the  outside  boundary,  r-r^  or  n“llmaxf 

Now  if  at  r=ro, 


Bm  dm 

3<a1,)  dJ™ 

and  so  on 

3r  dr 

and  if  we  force; 

H  cx)s(m^*))  =  0  H(A^(^’,'n)  Sin(m4‘))  =  0 

then  Es(^,T|)  matches  into  the  cylindrical  coordinate  system 
smoothly.  The  operator  H  is  the  Helmholtz  equation  in 

space . 

H(v)  =  a0-2P^  +  r0+J^I<^/  =  O 

The  functions  .  .  .  are  generalized  eigenfunctions,  which  are 

determined  numerically  by  solving  the  H  equation  subject  to 
proper  boundary  conditions.  At  T|*l  A^^,  .  .  .  are  simply 


constants . 


at  Tl-T|„ax 


dAj  raJJkr)^ 

dr  I  dr 

The  solution  for  A^„  consists  of  a  two  point  Boundary  value 
problem  which  requires  a  shooting  method  under  general 
conditions . 

As  an  illustration,  consider: 

H(A;i,  cos(m^‘))  =  0 


The  solution  procedure  involves; 

1.  Set  (aJ,  cos(m^*))^s^^  =  JJkTo)  cos(m^‘) 

2 .  Pick  an  initial  value  for  (A^)t^«1  =  constant  =  a; 

3 .  Solve  H  [A^  cos(m4 )]  =  0 

4 .  Check  the  slope  at  r=ro  to  see  whether  or  not : 

dr  dr 


If  agreement  is  unsatisfactory,  pick  a  new  value  of  A*^  and  go 
back  to  step  2 . 

In  turn,  all  four  functions  can  be  generated  and  four  general 
constants,  (A^n,)*,  (A^n,)*,  which  depend  only  on  the  problem 

geometry,  are  found.  Now,  we  can  match  the  B.C.  on  the  obstacle 
at  11=1  by  using  the  following: 
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■K-irAj, 
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Cm 

CM  E 
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+b2 

■K-irB^ 

1 

E 
>Q 
_ 1 

.-Ai, 

-(-I)”  Ail, 

-Bii, 

am 


a. 


-m 


Solving  easily  for  a^,  a_„,  b-„. 

Note  that  the  incident  angle,  (j>o,  affects  only  the  RHS  values. 
Therefore,  <|>o  can  be  varied  without  needing  to  recalculate  the 
eigenfunctions  .... 


Finally,  for  m=0, 

Aq  =  a^  Aq  —  bg  Bg 


Cg  =  ag  BJ  +  bg  aJ 


as  r  ro 


Aq  Jg(krg) 

aA;  i^a(Jgkr)^ 
ar  ar 

h(a;)=o  and  analogously  for 
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APPENDIX  D 


INFINITELY  LONG  CIRCULAR  CYLINDERS  -  NORMAL  INCIDENCE 
I.  Formulation 

Scattering  from  an  infinitely  long  cylinder,  when  the 
incident  field  propagates  normal  to  the  cylinder  axis  (taken  as 
the  z  axis)  is  essentially  a  two  dimensional  problem.  This  is 
evident  because  there  is  no  variation  in  the  incident  field  or 
in  the  cylinder  properties  in  the  axial  or  z  direction;  hence, 
the  scattered  field  will  not  vary  in  the  z  direction.  The 
cylinder  surface  is  a  natural  coordinate  surface  in  this  system 
(i.e.,  r  =  a),  and  boundary  value  techniques  give  an  exact 
solution . 

The  incident  wave  is  considered  to  be  planar,  and,  in 
general,  is  polarized  arbitrarily  r  jwever,  it  can  be  resolved 
into  two  components  :  (1)  an  E-field  component  parallel  to  the 
cylinder  axis,  and  (2)  an  E-field  perpendicular  to  the  same  axis 
(or  an  H-field  component  parallel  to  the  axis)  .  Each  case  is 
then  treated  separately.  The  two  cases  will  be  distinguished  by 
I  or  or  TM  and  TE,  referring  to  (1)  and  (2)  above 
respectively.  The  plane  of  incidence  is  considered  to  be  the  x-z 
plane . 

Since  an  electromagnetic  field  can  be  expressed  entirely  in 
terms  of  either  the  electric  or  the  magnetic  field,  the  problems 
can  be  solved  for 


Plane 

incident 

wave 

Fig.  1  Scattering  gaomatry  for  an  infinitely  long  circular  cylinder  with  axis  parpandicular  to  tha  page. 
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the  two  cases  by  :  (a)  using  only  the  axial  electric  field  for 
case  (1),  and  (b)  using  only  the  axial  magnetic  field  for  case 
(2)  .  In  either  case,  both  the  incident  and  scattered  fields 
posses  only  axial  components.  The  boundary  conditions  at  the 
cylinder  can  be  expressed  in  terms  of  these  axial  components, 
and  a  scalar  problem  rather  than  a  vector  problem  must  be 
solved,  which  simplifies  the  solution  considerably.  (Such  is  not 
the  case  with  a  sphere,  because  the  incident  field  strikes  the 
sphere  at  different  angles  corresponding  to  different  points  on 
the  sphere;  hence,  that  problem  retains  a  vector  nature.) 

Consider  the  incident  wave  to  be  resolvable  into  components 
of  the  form 


-ik^x 


(d.l) 


where  the  field  is  incident  from  the  left  in  figure  1.  The 
scattered  field  at  point  P'  can  be  expressed  for  each  of  the 
above  two  scalar  field  components  as 

oo 

e;(P0  =  Eji^^HfEnAnH^^koDcosCnf )  (d  •  2 ) 

oo 

Hl(P0  =  Hji.2^(-i)"£nBnH?Kr)cos(nf )  (d  •  3 ) 


and 

niO 

where 

hV(po= 

-lei 

and 

en=M 

for  n  =  0 

[2  for  n  =  1,2,3,... 

In  the  above  equations,  (^o^)  is  a  cylindrical  Hankel 

function  of  the  first  kind  of  order  n  and  argument  kor;  it  has 
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the  proper  behavior  large  r  corresponding  to  the  time  convention 
used  here.  The  constants  An  and  Bn  are  functions  only  of  the 
cylinder  material  and  radius,  and  are  given  in  following 
sections  for  the  cylinder  type  considered.  It  is  again  pointed 
out  that  in  case  (1)  ,  both  the  incident  and  scattered  electric 
fields  are  polarized  parallel  to  the  cylinder  axis;  while  in 
case  (2),  the  incident  and  scattered  magnetic  fields  are 
polarized  parallel  to  the  axis. 

II.  Far-Zone  Simplification  and  Scattering  Width 

Simplification  of  the  above  series  solution  is  possible 
when  (i)  kor  >>  (kor)^  and  (ii)  kor  >>  1  (or  r  »  Xq  )  •  this 

case,  a  large-argument  expansion  can  be  used  for  the  Hankel 
function.  That  is 

) 

where 

oe 

r>»C 


E.!(P0 


e 


T|,(f)  (d.4) 

) 

-T_L(f)  (d.5) 


and 


n-0 


In  contrast  with  the  exact  Mie  series  for  far-zone 
scattering  from  a  sphere,  the  field  from  an  infinite  cylinder 
decreases  as  1/Vr,  rather  than  the  1/r  dependence  for  objects  of 
finite  dimensions,  such  as  spheres.  Such  must  be  the  case  in 
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order  for  the  total  energy  in  the  diverging  scattered  fields  to 
be  conserved. 


( 

(d.7) 


and 


oi(*')  =  ^l  I* 

Two  special  scattering  directions  give  further  simplification  : 

(a)  Backscattering.  Here  ,  ({>'  =  0,  and  the  above  results 
simplify 

oo 

n=0 

and 

OO 

naO 

The  backscattering  widths  become 


=  I* 

^  n=0 

and 

oo 

<jl(0)  =  -i|  I* 

^  n-0 

(b)  Forward  Scattering.  Here,  =  7t,  and  the  scattering 
aunplitudes  simplify  to 


T„(7i)  = 


and 
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The  forward  scattering  widths  then  become 


af(7C)=-^|  1 2 

^  n=0 

and 


III.  Perfectly  Conducting  Cylinder 


The  coefficients  An  and  Bn  appearing  in  the  exact  solution 
for  the  perfectly  conducting  infinitely  long  cylinder  at  normal 
incidence  are 

A 

and 

g  Jn^M) 

H?^'(koa) 

where  the  primes  represent  derivatives  with  respect  to  the 
expressed  argument , Jn ( koa)  is  the  usual  cylindrical  Bessel 
function  of  order  n  and  argument  koa.  These  coefficients  are 
easily  generated  on  a  computer  by  recurrence  techniques,  or  are 
found  in  the  tables . 


